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Abstract 


This  report  presents  an  image  processing  technique  to  detect  satellite  streaks.  This  method  is 
particularly  useful  in  the  context  of  surveillance  of  space  where  the  positions  of  active 
satellites  and  other  orbital  debris  must  be  monitored.  In  such  cases,  the  orbital  parameters  are 
known,  but  after  a  certain  time  this  knowledge  deteriorates  because  of  the  recent  orbit 
perturbations.  Consequently,  the  satellites  need  to  be  re-observed  and  the  orbital  parameters 
must  be  updated  before  these  satellites  are  declared  lost. 

Hence,  even  before  its  reacquisition,  the  satellite’s  motion  is  known  but  its  position  is 
estimated  with  a  certain  margin  of  error.  This  knowledge  allows  for  automatic  pointing  of  the 
telescope,  acquiring  of  images  with  the  satellite  streak  within  the  field  of  view,  developing  a 
matched  filter  (with  the  motion  knowledge)  for  detecting  this  streak  within  the  image  and 
declaring  of  the  satellite  position. 

This  study  begins  with  a  detailed  analysis  of  a  typical  image,  which  includes  several  sensor 
artefacts  (such  as  dead  pixels,  background  gradient,  noise)  and  signal  degradation  (bleeding, 
blooming,  saturation,  etc.).  This  study  explains  how  sensor  artefacts  are  corrected,  image 
background  is  removed  and  noise  is  partially  erased.  Then,  it  describes  a  technique  to  detect 
and  erase  stars.  This  reduces  the  number  of  objects  in  the  image  that  could  generate  false 
alarms.  Finally,  this  document  describes  how  to  design  and  apply  a  matched  filter  used  for 
extracting  the  satellite  streak  in  images.  Examples  of  processed  data  are  illustrated  for  each  of 
the  processing  steps. 

Resume 

Ce  rapport  presente  une  technique  de  traitement  d’images  pour  detecter  les  traces  de  satellites. 
Cette  methode  est  particulierement  utile  dans  le  domaine  de  la  surveillance  de  Tespace  ou  les 
positions  des  satellites  actifs  et  d’autres  debris  doivent  etre  controlees.  Dans  ces  cas,  les 
parametres  orbitaux  sont  connus,  mais  apres  un  certain  temps  cette  connaissance  se  deteriore  a 
cause  des  recentes  perturbations  orbitales.  Consequemment,  les  satellites  doivent  etre  re¬ 
observes  et  les  parametres  orbitaux  mis  a  jour  avant  que  ces  satellites  ne  soient  declares 
perdus. 

Ainsi,  avant  meme  sa  re-acquisition,  le  deplacement  du  satellite  est  connu  mais  sa  position  est 
evalue  avec  une  certaine  marge  d’erreur.  Cette  connaissance  permet  de  pointer 
automatiquement  un  telescope,  d’acquerir  des  images  avec  la  trace  du  satellite  dans  le  champ 
de  vision,  de  developper  un  fibre  adapte  (avec  la  connaissance  du  deplacement),  de  detecter 
cette  trace  dans  T image  et  de  declarer  la  position  du  satellite. 

Cette  etude  commence  par  T analyse  detaillee  d’une  image  typique,  laquelle  inclut  plusieurs 
artefacts  de  capteur  (tel  que  des  pixels  morts,  gradient  de  fond,  bruit)  et  degradation  du  signal 
(coulage,  eblouissement,  saturation,  etc.).  Cette  etude  explique  comment  les  artefacts  du 
capteur  peuvent  etre  corriges,  le  fond  de  f  image  enleve  et  le  bruit  partiellement  efface.  Puis, 
elle  decrit  une  technique  pour  detecter  et  effacer  les  etoiles.  Cela  reduit  le  nombre  d’objets 
dans  Timage  qui  pourraient  generer  des  fausses  alarmes.  Finalement,  ce  document  decrit 
comment  concevoir  et  appliquer  le  fibre  adapte  qui  est  utilise  pour  extraire  la  trace  du  satellite 
de  Timage.  Des  exemples  de  donnees  traitees  sont  illustres  pour  chacune  des  etapes  du 
traitement. 
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Executive  summary 

This  report  describes  an  image  processing  technique  used  for  the  automatic  detection  and 
extraction  of  satellite  streaks  from  astronomical  images.  This  methodology  has  been 
developed  expressly  for  the  processing  of  images  produced  by  the  Surveillance  of  Space 
(SofS)  Concept  Demonstrator  (CD),  which  is  comprised  of  three  observatories  located  at 
Kingston,  Suffield  and  Valcartier.  Since  these  observatories  should  produce  hundreds  of 
images  every  week,  all  the  acquisition  and  data  reduction  tasks  must  be  automated  as 
much  as  possible. 

The  algorithms  presented  in  this  document  are  a  first  step  toward  a  fully  automated  data 
analysis  system.  They  automatically  process  and  remove  the  sensor’s  artefacts,  reduce 
the  image  content  by  removing  irrelevant  objects  and  finally  extract  the  satellite  streaks. 
This  last  step  is  possible  because  the  data  processing  is  integrated  with  the  tasking  and 
image  acquisition  system,  which  provides  a  priori  information  about  the  target  motion 
and  sensor  pointing.  It  only  lacks  a  detection  validation  step  to  the  system,  which  will  be 
the  object  of  future  work. 

These  image  processing  algorithms  were  implemented  within  the  operating  software  at 
the  Surveillance  of  Space  Operation  Center  (SSOC)  at  DRDC  Ottawa,  which  remotely 
operates  the  three  observatories.  Furthermore,  these  automatic  streak  extraction 
algorithms  should  also  be  implemented  into  the  SAPPHIRE  data  processing  facility  in 
the  near  future. 


Levesque  M.  P.,  Buteau  S.  2007.  Image  processing  technique  for  automatic  detection  of 
satellite  streaks.  DRDC  Valcartier  TR  2005-386.  Defence  R&D  Canada  -  Valcartier. 
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Sommaire 


Ce  rapport  decrit  une  technique  de  traitement  d’images  pour  la  detection  et  l’extraction 
automatique  des  traces  de  satellites  dans  les  images  astronomiques.  Cette  methodologie 
a  ete  specialement  developpee  pour  le  traitement  des  images  produites  par  le 
demonstrateur  de  concept  (CD)  de  la  surveillance  de  l’espace  (SofS),  lequel  est  constitue 
d’une  serie  de  trois  observatoires  localises  a  Kingston,  Suffield  et  Valcartier.  Puisque 
ces  observatoires  devraient  produire  des  centaines  d’images  chaque  semaine,  toutes  les 
taches  d’ acquisition  et  de  reduction  des  donnees  doivent  etre  automatisees  dans  la 
mesure  du  possible. 

Les  algorithmes  presentes  dans  ce  document  constituent  une  premiere  etape  vers  un 
systeme  entierement  automatise  d’analyse  des  donnees.  Ils  traitent  et  enlevent 
automatiquement  les  artefacts  des  capteurs,  reduisent  le  contenu  des  images  en 
supprimant  les  objets  non  pertinents  et  fmalement  ils  en  extraient  les  traces  de  satellites. 
Cette  demiere  etape  est  redue  possible  parce  que  le  traitement  des  donnees  est  integre 
avec  les  systemes  de  planification  des  taches  et  d’acquisition  d’images,  lesquels 
procurent  de  l’information  a  priori  sur  le  du  mouvement  de  la  cible  et  du  pointage  du 
capteur.  II  ne  manque  plus  au  systeme  qu’une  etape  de  validation  de  la  detection, 
laquelle  fera  l’objet  de  travaux  fiiturs. 

Ces  algorithmes  de  traitement  des  images  ont  ete  implantes  dans  le  logiciel  d’operation, 
au  centre  d’operation  de  la  surveillance  de  l’espace  a  Ottawa,  lequel  telecommande  les 
trois  observatoires.  De  plus,  ces  algorithmes  d’extraction  automatique  des  traces 
devraient  etre  implantes  dans  l’unite  de  traitement  des  donnees  de  SAPPHIRE  dans  un 
proche  avenir. 


Levesque  M.  P.,  Buteau  S.  2007.  Image  processing  technique  for  automatic  detection  of 
satellite  streaks.  DRDC  Valcartier  TR  2005-386.  Recherche  et  developpement  pour  la 
defense  Canada  -  Valcartier. 
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1.  Introduction 


In  the  current  state  of  colonization  of  near  Earth  space  by  satellites,  there  is  an  increasing 
need  to  know  exactly  the  real  status  of  occupation  of  this  space.  Thus,  orbital  parameters 
for  all  objects  traveling  in  this  space  must  be  known  with  a  high  degree  of  accuracy,  and 
this  knowledge  must  be  periodically  updated  because  this  situation  is  always  changing. 
Atmospheric  drag,  solar  wind,  moon  and  planetary  gravitational  perturbations,  Earth 
oblateness,  etc.  are  all  sources  of  interference  that  generate  orbital  perturbations  beyond 
what  the  best  orbital  model  can  predict.  The  solution  is  to  periodically  observe  all  the 
satellites,  particularly  the  debris  (because  active  satellites  themselves  contribute  to 
maintain  the  knowledge  of  their  orbital  parameters),  determine  with  precision  their 
positions  and  update  their  known  orbital  parameters. 

To  achieve  this  mission  of  Surveillance  Of  Space  (SOS),  the  Space  Surveillance 
Network  (SSN)  uses  several  large  telescope  and  radar  stations  located  all  around  the 
world.  However,  these  expensive  observation  stations  do  not  suffice  for  the  observation 
needs,  with  the  result  that  objects  are  lost  (due  to  not  being  reacquired  in  time)  every 
week.  At  the  same  time,  all  the  observations  do  not  really  require  execution  by  such 
powerful  sensors.  In  order  to  help  to  accomplish  the  mission,  a  battery  of  small  and  low- 
cost  instruments  (almost  amateur  grade  equipment)  can  be  deployed  to  contribute  to  the 
necessary  data  acquisition  tasks.  For  this  purpose,  the  US  has  developed  the  RAVEN 
stations  (Refs.  1  and  2)  and  Canada  also  contributes  to  the  SSN  with  similar  stations 
(Refs.  3  to  6),  which  were  initially  called  CASTOR  (Canadian  Automated  Small 
Telescope  for  Orbital  Research)  for  the  prototype  and  finally  rename  the  SofS/CD  for 
‘Surveillance  of  Space  Concept  Demonstrator’  for  the  operational  deployed  stations. 

A  CASTOR  station  (Refs.  3  to  8)  is  based  on  a  Celestron  14  telescope  coupled  with  a 
temperature  compensating  focuser,  mounted  on  a  Paramount  1 100GT  German- 
Equatorial  mount  (manufactured  by  Software  Bisque),  equipped  with  an  Apogee  AP8p 
1024x1024  pixel  CCD  camera.  This  instrument  is  located  inside  a  10’6”  dome 
manufactured  by  Ash  Domes  and  positioned  by  a  dome  control  system  manufactured  by 
Meridian  Controls.  Three  of  these  stations  are  spread  (to  reduce  the  risk  of  having  a 
cloud-covered  sky)  over  the  Canadian  territory  (at  Suffield,  Alberta,  Kingston,  Ontario 
and  Valcartier,  Quebec)  and  are  remotely  operated  from  a  common  operations  room 
(SOC:  Sensor  Operation  Control)  located  at  DRDC  Ottawa. 

When  the  certification  tests  are  completed,  these  three  stations  will  operate  in  automatic 
mode.  During  the  day,  the  operator  will  establish  the  observation  program  for  the 
following  night.  He  will  be  helped  by  a  schedule  optimization  program  (also  in 
development).  The  program  will  be  executed  at  night.  The  next  day,  the  operator  will 
check  system  errors,  reprogram  the  stations  and  communicate  the  valid  observation  to 
the  SSN.  This  implies  that  the  dome  opens,  rotates  and  closes  by  itself  (functions 
assumed  by  the  Meridian  dome  controller),  the  telescope  slews  to  the  appropriate 
position  at  the  right  time  (functions  that  can  be  assumed  by  the  Paramount  1 100GT 
German-Equatorial  mount  along  with  the  software  developed  by  Software  Bisque,  i.e., 
‘The  Sky’,  Automadome,  Orchestrate  and  CCD  Soft,  etc.)  and  that  the  camera  acquires 
the  images  at  the  right  time  (a  GPS  receiver  provides  an  accurate  clock).  The  whole 
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cycle  takes  less  than  3  minutes  for  one  satellite.  The  subsequent  (longer)  step  is  to 
download  the  image  from  the  CCD  camera  to  the  computer.  Hence,  during  good 
atmospheric  conditions,  it  is  foreseen  that  a  single  station  will  have  the  capability  to 
perform  up  to  300  observations  during  a  long  winter  night.  The  issue  lies  in  analyzing 
all  these  images. 

There  is  a  need  to  develop  algorithms  and  software  that  can  automatically  detect  and 
report  the  presence  of  satellite  streaks  in  the  acquired  images.  The  algorithms  presented 
in  this  document  were  developed  for  this  purpose.  They  were  tested  and  since  they 
performed  very  well,  they  were  incoiporated  into  the  remote-operation  software  at  the 
SOC.  They  are  also  efficient  regarding  processing  requirements.  The  first  satellite 
streak  detection  software  showed  that  the  unused  CPU  time  available  in  the  3  minute 
observation  cycle  was  sufficient  to  complete  the  analysis  task.  In  this  3  minute  cycle,  the 
computer  can  control  the  station,  analyze  the  last  acquired  image  and  report  on  the 
satellite  observation. 

The  image  processing  technique  presented  in  this  document  is  a  collection  of  algorithms 
used  to  detect  and  classify  everything  that  can  be  observed  in  the  image,  such  as  stars, 
satellite  streaks  and  image  artefacts.  The  basic  approach  consists  in  using  a  matched 
filter  to  detect  the  streak.  However,  since  it  cannot  be  detected  directly  with  high 
efficiency,  the  applied  technique  consists  in  detecting  anything  that  can  be  identified  as 
non-streak  (based  on  morphological  characteristics)  and  then  in  removing  it  from  the 
image,  leaving  a  simpler  image  where  the  streak  is  easier  to  detect.  This  process  reduces 
the  probability  of  false  alarms  and  increases  the  sensitivity  of  the  entire  algorithm 
system. 

This  report  begins  with  an  analysis  of  a  typical  image  where  signal  components  are 
characterized.  Next,  special  algorithms,  for  the  extraction  of  each  one  of  these 
components,  are  described.  One  by  one,  all  the  non-streak  components  (in  order:  sensor 
artefacts,  background,  noise  and  finally  the  stars)  are  eliminated  from  the  image  before 
the  streak  detection  algorithm  is  applied  to  the  remaining  signal.  All  previous  algorithms 
were  designed  to  preserve  streaks,  i.e.,  if  in  doubt,  the  processing  scheme  leaves  the 
signal  intact  to  avoid  altering  streaks,  even  though  the  output  signal  may  still  contain 
non-erased  stars. 

This  report  describes  several  inter-dependent  ad-hoc  algorithms.  These  algorithms  are 
inter-dependent  in  the  sense  that  some  of  them  are  specifically  designed  not  to  inhibit  the 
action  of  the  next  processing  steps.  In  other  words,  they  were  optimized  on  the  basis  of 
fine  trade-offs  between  the  different  processing  steps  rather  than  through  an  independent 
step-optimization  procedure  that  would  not  take  into  consideration  the  negative  impacts 
due  to  algorithm  interactions.  For  instance,  a  star  erasing  algorithm  with  a  lower 
threshold  would  be  more  efficient  for  star  removal  but  it  would  also  erase  faint  streaks. 

In  such  a  case,  it  is  preferable  to  tolerate  a  few  faint  stars  (which  are  tolerated  by  the  next 
processing  step)  and  preserve  the  faint  streaks.  Globally,  the  algorithms  presented  in  this 
report  are  individually  sub-optimal  and  consequently  could  certainly  be  further 
optimized  in  the  near  future.  Nevertheless,  it  is  shown  that  their  efficiency  and 
performance  were  adequate  to  the  current  application. 
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It  would  be  interesting  to  compare  the  performance  of  the  algorithms  developed  in  this 
report  with  other  known  methods.  Unfortunately,  there  are  almost  no  publications  on  this 
subject  in  the  open  literature.  This  is  probably  because  the  images  were  always 
examined  by  human  observers  and  the  conception  of  fully  automated  (operation  and  data 
processing)  is  only  a  recent  requirement.  One  may  think  that  detecting  artificial  satellites 
and  asteroids  (currently  done  by  astronomers)  are  similar  process,  but  unfortunately  the 
data  acquisition  constraints  are  different  and  the  characteristics  of  generated  data  are  not 
equivalent.  For  artificial  satellite  detection,  only  one  paper  deserved  to  be  mentioned. 
This  paper  by  Sanders-Reed  (Ref.  9)  describes  a  maximum  likelihood  filter  technique 
based  on  an  accurate  noise  model.  Although  the  maximum  likelihood  filter  is  certainly 
better  at  detecting  other  ‘unexpected’  satellites,  a  matched  filter  technique  was  favoured 
in  the  current  approach  because  of  its  higher  sensitivity  for  the  ‘expected’  satellites. 

The  detection  method  developed  in  this  report  is  implemented  in  the  SOC  and  has  been 
tested  with  hundreds  of  images.  It  has  proven  to  be  very  stable  and  sensitive.  The 
results  showed  that  it  is  able  to  detect  streaks  with  a  signal-to-noise  ratio  above  two.  The 
exact  figures  of  merit  are  not  evaluated  yet. 

At  the  moment  this  report  is  being  published,  the  detection  algorithms  are  being  refined 
(for  an  higher  sensitivity)  with  the  addition  of  a  post-detection  false-alarm  rejection 
algorithm  and  the  production  of  these  figures  of  merit  has  begun.  This  will  be  published 
in  Ref.  10. 

This  work  was  performed  between  September  2004  and  February  2005  under  the  project 
1 5 et  1 3 :  ‘Small  telescope  for  Surveillance  of  Space’. 
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2.  Image  characterization  and  correction 


The  first  step  in  the  development  of  an  image-processing  algorithm  is  to  analyze  typical 
images  in  order  to  characterize  all  signal  components.  For  this  purpose,  let  us  begin  with 
a  typical  case  such  as  the  one  below. 

Figure  1  shows  a  typical  image  acquired  with  the  CASTOR  system.  In  the  middle  of  the 
image,  there  is  a  very  small  and  faint  satellite  streak.  By  inspecting  the  digital  image 
along  with  contrast  manipulation  software,  a  trained  user  will  detect  it  without  difficulty. 
However,  the  system  is  designed  to  acquire  hundreds  of  images  per  clear  night  and 
fatigue  could  become  a  factor  that  would  limit  the  operator’s  performance.  Therefore, 
there  is  a  serious  need  for  an  algorithm  that  can  automate  the  detection  process,  i.e.,  read 
the  acquired  image  and  generate  a  formatted  report,  which  could  be  used  as  an  input  by 
another  program.  This  report  could  contain  the  coordinates  of  the  streak  end-points  along 
with  a  very  accurate  time  of  reference  and  the  astronomical  positions. 


Figure  1.  Typical  image,  acquired  with  amateur-grade  telescope  and  CCD,  that  contains  a 

faint  satellite  streak 

The  goal  of  the  processing  method  presented  below  is  to  extract  such  a  very  faint  signal 
almost  lost  in  the  noise,  among  sensor  artefacts  and  strong  signals  of  bright  stars.  It  is 
mandatory  to  carefully  analyze  images,  identify  signals  and  artefacts  and  make  a  good 
image  characterization.  Once  this  is  done,  one  can  note  that  these  image  components  are 
separable  and  can  be  removed.  After  removing  these  “extraneous”  image  components 
(bright  stars,  sensor  artefacts),  the  faint  satellite  streak,  almost  lost  in  the  original  image, 
becomes  more  readily  visible.  The  processing  begins  with  identifying  each  of  the 
potential  sources  of  trouble  and  with  finding  an  adequate  solution  to  each  one  of  them. 

In  the  following  sections,  all  the  image  artefacts,  which  are  taken  into  account  in  the 
whole  process,  are  presented.  Some  of  them  are  obvious  and  do  not  really  need  detailed 
descriptions.  Others  are  more  complex  and  require  special  attention. 
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2.1  Camera  artefacts 


The  first  noticeable  image  artefact  (Fig.  19-B1)  is  the  gradient  of  the  image  background. 
This  artefact  is  very  important  and  requires  special  attention  but  it  does  not  really 
represent  a  big  challenge.  This  subject  will  be  described  in  detail  later. 

The  second  important  image  artefacts  are  the  detector’s  bad  pixels.  This  problem  has  a 
very  classical  solution,  which  consists  in  interpolating  the  known  bad  pixels  using  the 
values  of  its  neighbours.  This  method  only  requires  building  a  list  of  bad  pixels  (or  line 
segments),  which  is  specific  to  each  CCD  camera  but  common  to  all  images  acquired  by 
the  same  CCD. 

Another  apparent  acquisition  artefact  is  the  dark  window  border  effect.  Very  often,  the 
last  image  lines  and  columns  (sometimes  up  to  20  consecutive  lines)  are  very  dark.  After 
analysis,  it  seems  that  the  gain  factor  is  incorrect  for  these  lines.  This  image  artefact  is 
compensated  by  using  a  polynomial  fit  on  the  image  background  (the  same  method 
described  below  for  the  background  removal  method)  and  by  checking  whether  the  last 
lines  (and  columns)  have  average  values  lower  than  expected.  For  these  image  lines,  the 
gain  is  readjusted  and  the  border  artefact  disappears.  The  elimination  of  this  dark 
window  effect  makes  the  signal  more  homogeneous  and  prevents  further  difficulties  for 
the  next  processing  steps. 


2.2  Image  noise 

Finally,  the  last  sensor  artefact  and  not  the  least  is  the  noise,  principally  generated  at  the 
acquisition  by  the  CCD  dark-current  noise  (almost  Gaussian),  the  readout  system  and 
partially  by  the  photon  noise  (Poisson  noise).  When  evaluating  the  image  background, 
the  CCD  noise  is  dominant  (at  least  for  the  amateur-grade  CCD  that  was  used).  This 
noise  directly  depends  on  the  sensor  temperature.  With  the  limited  performance  of  the 
cryo-cooler  install  on  the  camera,  very  often  the  CCD  is  not  cooled  enough. 

Furthermore,  the  noise  statistics  are  not  constant  in  the  image.  The  noise  is  at  its  lowest 
values  in  the  top  of  the  image  where  the  pixels  are  the  first  to  be  read  by  the  ADC 
(analog  to  digital  converter).  The  noise  level  gradually  increases  toward  the  bottom  of 
the  image  (Fig.  19-B2),  varying  from  the  simple  to  the  double.  This  indicates  that  there 
are  cumulative  effects  due  to  the  CCD  readout  mechanism.  Therefore,  the  local  noise  is 
always  a  factor  under  consideration  in  the  filter  design  presented  below.  It  is  used  as  the 
determining  criteria  to  establish  the  limit  of  sensitivity  for  the  local  filters.  For  this 
purpose,  the  standard  deviation  is  calculated  assuming  Gaussian  noise  (such  as  shown  in 
Fig.  2).  This  assumption  is  not  exact  but  it  is  accurate  enough  for  the  current  application. 

However,  from  time  to  time,  single  pixels  (or  pairs  of  pixels)  have  brightness  far  above 
the  Gaussian  distribution  (typically  10an  to  20an).  The  origin  of  these  noise  spikes  has 
not  yet  been  investigated.  They  are  not  necessarily  caused  by  hot  pixels  (bad  pixels  with 
an  excessively  high  gain)  because  they  do  not  occur  systematically  at  the  same  locations 
in  every  image.  They  may  simply  result  from  a  combination  of  several  noise  sources 
with  different  distributions.  It  is  certain  that  their  occurrence  has  a  very  negative 
influence  on  [local]  noise  estimation  and  on  the  star  detection  algorithms.  However, 
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they  are  very  easy  to  detect  and  eliminate.  When  detected,  they  are  set  to  zero  (or  at 
least  to  the  local  signal  average  if  the  image  background  has  not  yet  been  removed)  and 
the  remaining  image  background  is  easily  processed  as  if  the  remaining  noise  was  truly 
Gaussian. 


%JL  J 
Mil /I 

l.ll  1  1 ., 

Mill 

TIT 

1111  nmf 

|’  K||l  |  I'D  'Hf 

A:  Background  noise:  (almost 
Gaussian) 

Image:  signal  clipped  to  ±  2a  „. 


B:  Noise  spikes: 

Image:  signal  clipped  above  3a„. 

Spikes  are  single  pixel  with  exceptional 
brightness  far  above  main  distribution 
(»3an)  and  they  are  easily  identified  and 
removed. 


3an 

spike 

U/-A* 

Figure  2.  Gaussian-like  noise  and  noise  spikes 


2.3  Signal  dependent  artefacts 

The  next  category  of  artefacts  depends  on  the  sensor  response  to  the  signal  itself.  The 
signal  saturation  is  an  expected  phenomenon.  On  the  other  hand,  the  blooming  and 
bleeding  are  less  expected  and  must  be  carefully  examined  to  evaluate  their  impacts  on 
the  performance  of  the  streak  detection  algorithm.  Figure  3  illustrates  what  happens  to 
the  signal  properties  when  different  star  brightnesses  are  encountered. 

For  a  bright  star  (see  Figure  3)  that  is  well  above  the  noise  level  but  still  far  from  the 
saturation  level,  the  optical  impulse  response  is  clearly  visible  in  the  star  profile. 
Typically,  this  represents  a  star  that  covers  approximately  10%  of  the  sensor  dynamic 
range.  With  the  optical  setup  used  at  the  CASTOR  stations,  i.e.,  the  combination  of 
atmospheric  transmission,  telescope  and  CCD,  the  width  of  the  star  profile  at  half  height 
is  around  three  pixels.  This  is  the  star  radius  that  was  considered  in  the  filter  design. 
However,  for  other  systems  with  higher  resolution  (very  often  the  Nyquist  criterion  is  not 
applied  in  astronomical  systems)  this  value  should  be  different. 

For  fainter  stars,  the  star  size  remains  the  same,  except  that  the  profile  (the  PSF,  i.e.,  the 
Point  Spread  Function)  becomes  more  irregular  due  to  the  sampling  effect.  For  very 
faint  stars,  the  profile  cannot  be  used  but  the  aggregation  of  pixels  remains  visible  for 
object  brightnesses  only  a  few  sigmas  (standard  deviation  of  the  noise)  above  the 


6 


DRDC  Valcartier  TR  2005-386 


background.  These  last  objects  can  be  ignored  because  they  do  not  represent  a  limiting 
factor  for  the  satellite  streak  detection. 

On  the  opposite  end,  the  very  bright  stars  (see  Figure  3)  can  be  a  serious  source  of 
problems.  First,  for  the  streak  detection  algorithm,  they  are  sources  of  signal  that  are  in 
competition  with  the  satellite  streak.  Therefore  they  need  to  be  detected  first  and 
removed  from  the  image.  Second,  when  the  brightness  increases,  the  width  of  the  star 
profile  also  increases.  This  can  only  be  explained  by  a  sensor  deficiency;  when  heavily 
exposed,  the  electron  charges  in  the  CCD  pixels  increase  to  the  point  that  the  electrons 
are  repulsed  in  the  neighbour  pixels,  producing  a  blooming  effect.  Thus,  the  filter 
designed  to  detect  and  remove  the  stars  requires  several  passes,  where  each  pass  is 
adapted  for  a  specific  range  of  star  brightness  and  size. 

Afterwards,  when  the  signal  intensity  continues  to  increase,  the  problem  of  saturation 
appears.  The  streak-processing  algorithm  can  tolerate  a  few  slightly  saturated  stars. 
Flowever,  a  single  star  with  several  saturated  pixels  is  an  indication  of  other  fatal 
artefacts.  First,  the  half-height  width  of  the  star  profile  does  not  increase  proportionally 
with  the  number  of  saturated  pixels.  Flence,  it  is  important  to  count  these  saturated 
pixels.  A  count  above  a  certain  limit  is  an  indication  that  the  image  is  degraded  by 
severe  artefacts  (visible  diffraction  and  reflection  patterns)  and  this  image  should  not  be 
process.  It  should  be  simply  discarded.  Second,  the  aspect  ratio  of  these  saturated  pixels 
is  an  indication  of  another  problem:  bleeding.  Figure  3  shows  a  saturated  star  where  the 
bleeding  has  just  begun  to  appear  (aspect  ratio  smaller  than  2)  and  a  very  saturated  star 
with  a  catastrophic  bleeding.  Algorithms  were  developed  to  erase  the  bleeding  by 
interpolating  the  saturated  pixels  from  the  nearest  valid  pixels.  Flowever,  additional 
analyses  have  shown  that  this  interpolation  is  not  useful  when  severe  bleeding  is  detected 
(i.e.,  with  a  height/width  aspect  ratio  greater  than  2  or  3)  because  the  image  is  useless. 

As  a  matter  of  fact,  a  strong  bleeding  is  an  indicator  that  the  star  is  so  bright  that  several 
diffraction  patterns  appear  in  the  image  (halos,  secondary  and  parasitic  reflections,  etc.). 
For  that  case,  the  image  background  cannot  be  removed  (because  of  complex  structure  of 
the  star  pattern)  and  the  satellite  streak  cannot  be  detected.  Therefore,  it  is  important  to 
program  the  telescope  to  avoid  acquiring  satellite  streak  in  the  presence  of  very  bright 
stars.  Prior  to  processing  the  image,  saturated  stars  must  be  counted  and  if  the  sum 
exceeds  a  determined  threshold,  the  image  must  be  rejected. 
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Very  faint  stars: 

3an  <  SNR  <  10  an 

or  typically  signal  <  0.5%  of  the  sensor  dynamic  range. 

Comment:  Irregular  shape,  peak  profile  visible  on  5  to  9  adjacent  pixels, 
effects  of  the  optical  impulse  response  still  visible. 


Faint  stars: 

10an  <  SNR  or  typically  0.5%  to  4%  of  the  sensor  dynamic  range. 

Comment:  Gaussian  profile  obvious  due  to  the  optical  impulse  response 
(along  with  the  atmospheric  propagation),  star  completely  contained  in  a 
(6x6)  window 

(half-height  width  =  3  pixels) 


□ 

□ 


Medium  and  bright  stars: 

4%  to  15%  of  the  sensor  dynamic  range. 

Comment:  Clean  profile,  slightly  larger  impulse  response  caused  by  a 
light  blooming  effect. 

window  size  =  (8x8)  pixels  to  (10x10)  pixels 


Very  bright  stars: 

15%  to  100%  of  the  sensor  dynamic  range. 

Comment:  Strong  blooming  effect  visible,  requires  a  larger  window  size 
to  encircle  the  star:  (13x13)  pixels  and  more.  Furthermore,  the  Gaussian 
shape  of  the  star  has  a  long  decay  and  its  signal  is  perceptible  above  the 
background  noise  as  far  as  20  pixels  from  its  center,  even  if  its  half- 
height  width  is  only  7  or  8  pixels _ 


Saturated  stars: 

Comment:  Case  by  case,  local  window  size  depends  on  the  size  of  the 
saturated  area: 

-  4  saturated  pixels:  strong  blooming,  require  at  least  a  14X14 
window 

-  7  saturated  pixels:  very  strong  blooming,  requires  at  least  a 
25x25  window 


Very  Saturated  stars: 

Comment:  Gaussian  profile  severely  truncated  at  the  maximum  range. 
Very  large  area  of  saturated  pixels.  Not  only  the  blooming  artefact  is 
severe  but  there  are  also  strong  bleeding  artefacts.  When  detected,  there 
is  a  big  chance  that  the  image  cannot  be  used  to  detect  faint  satellite 
streaks  because  of  the  other  diffraction  artefacts  (halo)  introduced  in  the 
image.  The  image  shows  Vega  with  a  10  s  exposure  with  a  CCD 
mounted  on  a  40  cm  Celestron  telescope. 


Figure  3.  Star  morphology 
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2.4  The  satellite  streaks 

After  having  described  everything  that  we  want  to  get  rid  off,  it  remains  only  the  signal 
we  are  looking  for:  the  satellite  streak.  This  is  the  approach  used  to  facilitate  the  satellite 
detection.  It  is  difficult  to  detect  streaks  but  it  is  quite  easy  to  detect  everything  else. 
Hence  let  us  eliminate  everything  in  the  image  that  does  not  look  like  a  satellite  streak 
until  the  streak  becomes  one  of  the  dominant  remaining  objects. 

The  following  figure  illustrates  examples  of  satellite  streaks  that  the  algorithm  must 
detect.  If  streaks  were  always  very  bright,  a  simple  intensity  threshold  would  be  enough 
to  detect  them,  but  this  rarely  happens.  Since  it  is  desired  to  detect  fainter  and  fainter 
objects,  the  performance  of  the  detection  technique  must  be  optimized  to  detect  objects 
close  to  the  limit  imposed  by  the  noise  and  sensitivity  of  the  observation  system.  The 
algorithm  presented  in  this  document  can  detect  satellite  streaks  with  average  intensities 
lower  than  the  noise  standard  deviation. 

Usually,  the  satellite  streak  is  a  small  line  segment,  very 
often  almost  lost  in  the  noise.  Generally,  the  user  must 
manipulate  the  image  contrast  to  seen  the  satellite  streak. 
In  this  small  image,  the  intensity  of  this  streak  is  three 
times  the  noise  standard  deviation  and  less  than  the 
background  gradient  variation  (from  end  to  end). 


In  this  example,  the  average  intensity  of  the  streak  is 
0.75  a,  almost  lost  in  the  noise.  Due  to  the  great  number 
of  pixels  in  the  streak,  it  is  still  observable  and  detectable. 
However,  the  noise  level  prevents  determining  the  end¬ 
points. 


This  last  case  illustrates  a  very  clear  streak  (that  rarely 
happens).  It  also  illustrates  the  fading  effect  that 
frequently  occurs  with  low  orbit  objects,  because  of  their 
entry  into  the  Earth  limb.  Geostationary  and  Molnyia 
active  objects  are  usually  more  constant  (for  short 
exposure),  except  for  the  debris  that  often  rotates  and 
changes  its  reflective  properties  during  rotation. 


Figure  4.  Examples  of  satellite  streaks 
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3.  Processing  overview 


Now  that  the  image  is  properly  described,  it  is  easier  to  see  that  these  different  objects 
(noise,  background,  stars  and  streak)  are  separable  objects  that  can  be  processed 
individually.  Figure  5  illustrates  the  global  processing.  The  basic  idea  for  detecting 
difficult  streak,  as  already  described  previously,  is  to  first  remove  everything  else  that 
can  be  detected  and  removed. 

To  detect  an  object  with  a  known  shape,  the  maximum  sensitivity  is  achieved  with  a 
matched  filter.  However,  the  drawback  for  such  a  filter  is  the  generation  of  false  alarms 
by  other  objects.  Particularly,  the  stars  (almost  Dirac’s  delta)  produce  strong  matched- 
filter  responses.  Nevertheless,  the  stars  are  easy  to  detect  with  filters  that  use  local 
statistics,  but  such  filters  are  easily  eluded  by  the  image  background.  Fortunately,  the 
background  can  be  easily  evaluated  and  subtracted  using  polynomial  fits.  Therefore,  the 
background  is  removed  first,  followed  by  the  stars  and  finally  the  streak  is  detected,  as 
illustrated  in  Figure  5. 

Each  of  these  steps  requires  special  attention.  When  removing  the  background,  small 
features  above  the  average  background  (stars,  streaks,  etc.)  must  be  ignored.  Thus,  when 
a  pixel  signal  is  too  high,  comparatively  to  its  neighbours  (presence  of  a  star),  this  pixel 
is  not  considered  in  the  polynomial  fit  used  to  estimate  the  background.  This  is 
described  in  Chapter  5.  After  that,  before  the  stars  are  detected  using  local  statistics, 
special  attention  is  given  to  the  noise  spikes  (in  Chapter  6),  which  are  shaip  and  intense 
noise  events.  A  spike  is  electronic  noise,  and  since  it  is  not  created  by  photons 
propagated  through  the  optics,  it  does  not  have  the  typical  optical  impulse  response 
profile  (PSF).  It  is  narrower,  typically  only  one  pixel  width.  These  noise  spikes  severely 
reduce  the  performance  of  the  star  detection  algorithm.  Afterwards,  the  stars  are  easily 
detected  but  their  removal  is  tricky.  Their  profiles  must  be  evaluated  and  subtracted. 

The  details  are  explained  in  Chapter  7.  Finally,  the  matched  filter  is  applied  to  detect  the 
satellite  streak.  The  matched  filter  response  is  not  as  sensitive  as  desired,  but  an  iterative 
approach  can  meet  the  desired  detectability  levels.  This  method  is  presented  in  Chapter 
8. 
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Figure  5.  Processing  overview.  Images  at  various  processing  steps  are  shown  in  Fig.  1 9. 
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Detection  of  saturated  stars  and  bleeding 
estimation 


This  task  is  essential  to  evaluate  the  image  quality  because  when  extreme  bleeding  is 
detected,  the  image  is  probably  useless.  Therefore,  the  saturated  stars  need  to  be 
detected  and  characterized. 

The  process  begins  by  counting  the  total  number  of  saturated  pixels  (an  intensity  of 
65535  for  a  16-bit  image)  in  the  entire  image.  Then,  starting  with  a  saturated  pixel,  a 
region-growing  algorithm  is  applied  to  gather  all  the  contiguous  saturated  pixels  of  the 
first  star.  The  minimum  and  maximum  coordinates  in  height  and  width  are  recorded  and 
the  aspect  ratio  is  estimated.  After  that,  the  number  of  saturated  pixels  of  the  first  star  is 
subtracted  from  the  total  number  of  saturated  pixels  and  the  algorithm  searches  for  the 
next  star  until  there  are  no  more  saturated  pixels  in  the  image.  The  final  output  of  this 
algorithm  is  a  list  of  saturated  stars,  with  their  position,  number  of  saturated  pixels  per 
star  and  aspect  ratio.  This  list  will  also  be  used  later  for  the  first  estimation  the  image 
gradient  background,  because  it  indicates  non-background  pixels. 

This  approach  allows  differentiating  two  cases.  The  saturated  pixels  can  be  distributed 
into  a  few  severely  saturated  stars  or  into  several  “barely”  saturated  stars.  In  both  cases, 
the  total  number  of  saturated  pixels  can  be  the  same.  In  the  first  case,  each  star  has 
several  saturated  pixels  and  presents  bleeding  symptoms.  Images  of  this  kind  should  not 
be  processed;  they  are  useless.  For  the  second  case,  there  are  several  saturated  stars, 
each  one  having  only  a  few  saturated  pixels.  These  images  can  be  processed  correctly. 
Experiments  have  shown  that  the  streak  detection  algorithms  performs  well  with  several 
stars  containing  no  more  than  5  to  1 0  saturated  pixels  each.  However,  a  single  star  that 
has  more  that  25  saturated  pixels  (and  an  aspect  ratio  height/width  above  3  indicates  a 
bleeding  effect)  is  an  indication  that  the  CCD  is  bloomed  and  that  optical  artefacts 
(diffraction,  reflection,  etc.)  will  be  visible  in  the  image. 
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5.  Background  removal 


Background  removal  is  a  critical  step  in  the  current  processing.  The  performance  of  the 
star  and  streak  detection  algorithms  depends  on  the  quality  of  the  background-removal 
algorithm.  The  later  must  not  leave  any  residual  background  that  could  corrupt  the 
evaluation  of  local  noise  and  signal  statistics.  The  local  statistics  evaluation  is  based  on 
the  assumption  that  the  background  is  null. 

The  dark- frame  method  was  first  tested.  However,  two  consecutive  expositions  (with 
identical  condition,  even  two  dark  frames)  do  not  really  generate  identical  backgrounds. 
The  image  gradient  is  slightly  different  from  one  exposition  to  another  and,  when 
subtracted,  they  never  cancel  each  other  perfectly.  Experience  showed  that  this  is 
principally  due  to  cryo-cooler  deficiencies.  The  problem  could  have  been  avoided  by 
using  better  hardware  operated  in  better  conditions.  It  was  however  decided  to  operate 
the  cheapest  possible  hardware  (amateur  grade)  as  long  as  the  sensitivity  requirement  is 
met.  Moreover,  the  dark  frame  method  doubles  the  sensor  workload,  which  is  not 
desirable. 

Several  commercial  software  products  were  also  tested  to  remove  the  image  background 
gradient.  Usually,  these  software  packages  provide  significant  image  improvements  but 
some  image  artefacts  (residual  local  background)  remain  that  are  not  desirable  for  the 
streak  detection  purpose.  Therefore,  it  was  decided  to  develop  a  specific  method  for 
background  removal.  This  method  takes  into  account  the  characteristics  of  the  sensors 
used  in  the  system  and  the  application  constraints. 

This  method  is  based  on  the  fact  that  the  image  background  is  very  smooth  and  can  be 
modeled  with  polynomial  fits.  However,  the  presence  of  signals  (stars  and  streak) 
prevents  a  classical  polynomial  method  from  obtaining  exact  background  estimation  and 
the  result  of  a  blind  fit  is  inadequate.  Thus  the  method  is  improved  with  a  few  intelligent 
rules,  which  consists  of  an  initial  estimate  for  detections  of  stars  that  will  be  ignored  in 
the  background  estimation.  However,  to  detect  a  star  with  its  brightness  and  an  intensity 
threshold,  the  algorithm  needs  to  evaluate  the  background,  which  does  not  seems  make 
sense  since  the  problem  is  in  the  answer.  This  problem  is  solved  with  an  iterative 
process  such  as  shown  in  Fig.  6.  Therefore,  a  first  blind  iteration  is  done  where  the 
polynomial  is  fit  the  image  which  includes  the  stars.  Using  this  first  preliminary 
background  estimation,  a  second  polynomial  fit  is  calculated  where  all  objects  brighter 
than  the  previously  estimated  background  are  truncated  in  intensity  (using  a  rough 
threshold),  producing  a  better  estimation.  Finally,  this  same  operation  is  repeated  a  third 
time  with  a  threshold  that  is  well  adjusted  to  the  background.  This  third  estimation 
truncates  every  star  to  produce  estimation  satisfactory  for  our  purposes. 
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Figure  6.  Iterative  polynomial  fit  process  to  remove  the  image  gradient  background 


Before  describing  the  method  in  detail,  it  must  be  mentioned  that  the  image  columns  are 
easier  to  model  than  the  lines  (specifically  for  that  CCD  camera),  as  shown  by  the 
extracted  profiles  illustrated  in  Fig.  7.  These  profiles  illustrate  the  image  background 
gradient  along  with  star  inclusions.  One  can  see  that  the  vertical  gradient  is  basically  a 
‘ramp’  function,  which  is  probably  an  artefact  created  by  an  accumulation  of  bias 
electrons  in  the  analog  to  digital  converter  (ADC)  in  the  CCD  shift  register  mechanism. 
One  can  also  see  that  a  first  order  polynomial  can  easily  fit  that  background.  Moreover, 
a  careful  analysis  shows  a  small  third  order  component  that  should  be  taken  into  account. 
The  line  profiles  appear  to  be  second-order  curves  with  potentially  higher  order 
components  that  appear  harder  to  model.  Therefore,  it  was  decided  to  model  the 
background  with  polynomials  fitted  on  the  column  axis. 
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Figure  7.  The  column  and  line  profiles  extracted  from  the  image  indicate  that  the  background  is  by  far 

easier  to  model  over  the  column  axis 


As  written  above,  the  first  iteration  is  performed  without  clipping  the  stars,  except  for  the 
already  known  saturated  stars  (as  detected  in  the  previous  section).  The  areas 
surrounding  these  saturated  stars  are  ignored  in  the  ‘first’  background  estimation.  The 
sizes  of  these  areas  are  discussed  below  in  Section  7  ‘Star  removal’.  A  series  of  first 
polynomial  fits  of  the  third  order  are  then  calculated  over  the  ‘jth’  columns,  resulting  in  a 
rough  ‘first’  background  estimation  (pixel  (i,  j)  =  Poj  +  P|,  i  +  P2j  i2  +  P?,  i3 ).  After 
correction  with  this  ‘first’  estimation  (explained  in  next  paragraph),  a  ‘second’  improved 
background  estimation  is  obtained  by  ignoring  all  bright  stars.  This  time,  the  bright  stars 
are  simply  truncated  by  applying  a  coarse  threshold  value,  which  is  equal  to  the  previous 
background  estimation  plus  2%  to  5%  of  the  sensor  dynamic  range.  This  produces  a 
‘second’  background  estimation  with  an  error  of  only  a  few  counts  per  pixel.  This 
process  is  repeated  again  a  third  time  with  a  finer  threshold.  The  noise  level  (an)  is 
evaluated  and  the  fine  threshold  is  set  to  3  an  above  the  ‘second’  estimated  background. 
This  time  we  almost  obtain  a  perfect  ‘third’  background  estimation. 

Faint  horizontal  bright  bands  may  remain  at  this  step.  They  are  created  by  the  inability 
of  a  third  order  polynomial  to  exactly  model  the  background.  Using  a  higher  order  for 
the  polynomial  only  causes  a  change  in  the  pattern  of  these  horizontal  bands.  The 
polynomial  fit  generates  a  ‘bouncing’  effect  in  the  image.  This  artefact  is  eliminated  by 
using  the  same  polynomial  fit  again  over  the  line  axis  rather  than  the  column  axis.  Also, 
one  must  note  that  the  last  pixels  in  lines  or  columns  are  never  considered  in  the 
polynomial  fit  method  because  of  the  dark  edge  frame  artefact  described  above. 

As  mentioned  previously,  the  first  polynomial  fit  requires  special  attention.  This  first 
background  estimation  contains  ill-conditioned  column  fits,  due  to  the  presence  of  very 
bright  stars  that  strongly  influence  the  results.  If  they  are  not  taken  into  account 
immediately,  the  convergence  of  the  method  is  severely  reduced  and  may  require  several 
additional  iterations  before  the  background  is  correctly  removed.  However,  there  exists 
a  quick  way  to  immediately  locate  columns  that  contributed  to  bad  fit  results.  The 
polynomial  coefficients  of  every  column  polynomial  fit  are  compared  and,  since  the 
background  is  relatively  smooth,  coefficients  of  neighbouring  columns  should  be  very 
similar.  If  not,  the  differences  indicate  that  some  fits  are  corrupted,  as  shown  in  Fig.  8. 
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In  those  cases,  the  wrong  coefficients  are  replaced  by  values  interpolated  from  their 
neighbours.  These  ‘outlaw’  values  (in  blue  in  Figure  8)  are  simply  detected  with  a  low 
order  polynomial  fit  calculated  over  all  the  equivalent  coefficients  (e.g.  over  the 
coefficients  Po  of  every  fit  performed  for  every  image  column).  Then,  the  standard 
deviation  between  the  coefficients  and  the  result  of  that  fit  is  evaluated.  Finally,  the 
coefficients  whose  values  differ  by  more  than  3a,  relatively  from  the  main  sequence 
(indicated  in  blue  in  Fig.  8),  are  replaced  by  interpolated  values.  Lower  differences  (in 
red  in  Figure  8)  are  not  compensated  because  they  likely  represent  non-homogeneity  in 
the  CCD  channel  reading  mechanism  (each  image  column  represent  a  different  channel), 
and  it  is  a  desired  effect  that  the  background  removal  algorithm  also  removes  this 
artefact.  This  process  of  coefficient  comparison  needs  to  be  done  only  for  the  first 
iteration.  Afterwards,  the  background  removal  algorithm  converges  very  rapidly. 


Figure  8.  This  graph  shows  the  results  for  coefficient  P0  (but  the  graphs  are  similar  also  for  Pi, 
P2  and  P4)  after  the  first  iteration.  The  polynomials  were  calculated  over  each  of  the  1 024 
image  columns.  Since  the  background  is  very  smooth,  the  coefficients  should  be 
approximately  the  same  for  all  image  columns.  The  blue  spikes  indicate  the  columns  where 
incorrect  coefficients  Po  were  calculated  because  these  image  columns  contain  bright  stars. 
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6.  Noise  spike  removal 


Given  that  the  algorithm  for  star  detection  and  removal  is  based  on  a  local  estimation  of 
the  noise  and  signal,  it  is  very  important  that  the  noise  itself  does  not  include  erratic 
spikes  that  could  be  misinterpreted  as  a  real  signal.  During  the  execution  of  the 
background  removal  algorithm  (Figure  19B1),  the  noise  standard  deviation  ‘an’  (Figure 
19B2)  had  been  calculated  everywhere  in  the  image,  except  for  the  areas  containing 
stars.  At  these  locations,  the  noise  characteristics  were  interpolated.  Typically,  noise 
level  increases  by  (almost)  a  factor  two  from  the  top  to  the  bottom  of  the  image  for  the 
considered  CCD  camera.  Subsequent  noise  measurements,  with  a  small  window  (3x3), 
produce  local  estimations  ay,  which  are  roughly  equal  to  0.5  to  3  times  the  known 
standard  deviation  a„.  This  variation  is  normal  considering  the  limited  number  of  pixels 
(9  pixels)  for  the  estimation.  Flowever,  a  noise  spike  may  suddenly  produce  stronger 
variations,  sometimes  more  than  1 0  an. 

All  real  objects,  except  the  noise  spike,  have  a  width  imposed  by  the  optical  impulse 
response  (PSF:  point  spread  function).  As  already  shown  previously  in  Figure  3,  even 
the  faintest  stars  are  a  few  pixels  wide.  Figure  2B  shows  that  this  is  not  true  for  noise 
spikes  which  are  created  by  the  electronics.  Flence,  a  simple  but  efficient  filter  was 
designed  to  detect  and  remove  these  bad  pixels  (spikes)  that  could  generate  problems 
later. 

The  filter  principle  consists  in  erasing  the  intense  objects  whose  profiles  are  narrower 
than  the  optical  PSF.  For  a  PSF  with  a  half-height  of  3  pixels,  every  pixel  around  the 
pixel  of  maximum  intensity  (of  a  given  star)  has  at  least  half  the  intensity  of  this  pixel. 
For  a  noise  spike,  this  relation  is  not  preserved.  There  is  no  optical  convolution  that 
creates  dependence  between  the  pixels. 

The  filter  is  shown  in  Fig.  9.  When  a  pixel  ‘A’  shows  an  intense  signal  (>  3an)  and  its 
neighbour  pixels  (‘B’,  ‘C’,  ‘D’  or  ‘E’)  look  like  noisy  background  (<  2an),  then  the  pixel 
is  reset  to  zero.  This  removes  90%  of  the  noise  spikes. 

Flowever,  signal  analysis  showed  that  there  are  often  two  or  three  consecutive  noise 
spikes.  A  real  star  would  generate  a  two  dimension  array  of  pixels.  The  very  faint  star 
presented  in  Figure  3  has  approximately  5x5  pixels  and  the  PSF  function  is  visible  on 
both  axes.  The  image  contains  several  intense  objects  which  are  a  few  pixels  long  on 
one  axis  and  only  one  pixel  wide  on  the  other  axis.  Such  an  object  cannot  be  created  by 
photons  propagated  through  the  optical  system. 
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Single  pixel  spike: 

If  A  >  3a  and  (Bj<2a  and  Cj<2a  and  Dj<2a  and  Ei<2a)  :  A  =  0. 
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Two-pixel  spike: 

If  A  >  3a 

and  ( (B  i  >2a  and  B2<2a  )  xor 

( (Ci  >2a  and  C2<2a  )  xor  ...  up  to  E 
then  A=  0  and  B  or  C  or  D  or  E  =  0. 


Figure  9.  Geometry  of  the  kernel  used  to  remove  the  noise  spikes  like  those  illustrated  in  Fig.  2 


Therefore,  a  better  filter  can  be  conceived  by  considering  the  occurrence  of  consecutive 
spikes.  The  second  filter,  presented  in  Figure  9,  considers  that  an  intense  pixel  can  be  a 
spike  if  no  more  than  one  other  adjacent  pixel  is  more  intense  than  the  noisy  background. 
When  this  occurs,  both  pixels  are  reset  to  zero.  This  filter  removes  more  than  99%  of  the 
spikes  and  considerably  improves  the  performance  of  the  subsequent  filtering  operations. 
Considering  other  geometries  of  spike  occurrence  would  make  for  unnecessarily 
complex  algorithms  with  little  or  no  significant  performance  improvements. 
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7.  Star  removal 


The  streak  detection  algorithm  is  based  on  a  matched  filter  technique.  This  filter  readily 
detects  line  segment  shapes.  Its  main  drawback  is  the  false  alarms  that  might  be 
generated  by  other  objects  (like  stars).  Moreover,  bright  stars  (which  have  more  energy, 
or  counts,  in  their  peaks  than  a  satellite  streak  has  in  its  whole  trace)  may  generate 
stronger  responses  than  the  sought  satellite  streak.  Therefore,  a  straightforward  solution 
is  to  detect  and  erase  the  stars  prior  to  make  an  attempt  at  detecting  satellites. 

The  stars  are  easy  to  detect.  The  challenge  is  to  detect  and  erase  them  without  affecting 
a  yet  undetected  satellite  trace.  The  main  difficulty  lies  in  the  fact  that  different  stars 
have  different  characteristics,  as  shown  in  Fig.  3.  No  single  methods  can  erase  all  stars 
in  a  single  pass.  Therefore,  depending  on  the  star  brightness,  blooming  and  spatial 
distribution,  several  filtering  iterations  are  performed  with  filter  parameters  adapted  to  a 
specific  class  of  star. 


7.1  Star  detection: 

Stars  are  detected  and  erased  in  a  3-step  process,  starting  with  saturated  stars,  then  with 
bright  stars  and  finally  with  faint  stars.  The  saturated  stars  are  already  known  since  they 
were  detected  before  the  background  removal  process.  They  are  erased  using  the  method 
described  in  the  next  section. 

The  very  bright  and  faint  stars  are  both  detected  with  adapted  double-gate  type  filters. 
The  local-window  geometry  of  those  filter  are  adapted  to  the  star  apparent  radius 
(depending  on  the  blooming).  This  is  why  two  (if  not  three)  filter  passes  are  required 
with  different  window  geometry  and  threshold  values.  We  will  describe  in  the  following 
paragraphs  how  the  detection  filter  is  designed. 

Discriminating  stars  from  streaks  is  simple.  All  the  star  pixels  are  located  inside  a  very 
small  region,  whereas  the  streak  pixels  (the  direction  of  which  has  not  been  not  specified 
yet)  cannot  be  confined  in  that  same  region.  Hence,  using  a  double-gate  filter,  it  is 
possible  to  measure  the  star  signal  in  the  inner  window  (using  the  evaluated  values  of 
local  mean  p;n,  local  standard  deviation  am  and  the  local  maximum  intensity  maxm) 

while  the  outer  window  (with  pout  and  aout )  confirms  that  there  is  only  “background 
noise”  surrounding  the  star.  Figure  10  shows  that  this  is  true  only  for  a  single  star.  This 
is  false  for  a  streak  of  a  significant  intensity  because  there  is  also  a  measurable  signal  in 
the  outer  window. 

This  classical  double-gate  filter  is  convenient  when  the  star  profile  is  smaller  than  the 
inner  window,  but  this  is  not  always  the  case.  Furthermore,  the  outer  window  lacks 
sensitivity  because  when  it  tries  to  check  for  the  presence  of  the  streak  prolongation  in 
the  outer  window,  the  local  statistics  are  calculated  in  the  presence  of  a  large  area  of 
background,  which  reduces  the  filter  sensitivity.  Therefore,  this  basic  filter  is  modified 
and  its  new  design,  presented  in  Fig.  11,  proves  to  be  more  efficient  for  streak  detection 
in  the  outer  window. 
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Figure  10.  A  double-gate  filter  for  star  detection  where  the  9  pixels  in  the  inner  window  are  used  to 
evaluate  the  star  statistics  and  the  1 6  outer  pixels  evaluate  the  background 


Figure  1 1.  Multiple-window  filter  for  the  detection  of  faint  stars.  Inner  window  is  7x1  pixels  for  faint 
stars  and  could  be  increases  up  to  13x13  pixels  for  bright  stars ,  to  correct  for  blooming  effect.  The 
outer  windows  are  used  to  detect  the  presence  of  a  streak 

In  this  new  design,  a  faint  star  (with  a  half-height  width  of  3  pixels)  is  completely 
confined  inside  the  inner  7x7  pixel  window.  Its  average  value  pin  is  easy  to  evaluate 
over  this  area.  For  a  more  sensitive  method,  the  maximum-intensity  pixel  could  be 
found  (inside  that  area)  and  pm  could  be  evaluated  within  a  3x3  window  around  that 
point.  Hence,  only  the  brightest  pixels  would  be  considered.  Then,  in  the  attempt  to 
check  the  surrounding  background  (to  see  if  the  extended  object  can  be  a  streak  rather  a 
star),  a  new  configuration  including  several  small  outer  windows  (12  outer  windows  in 
the  case  of  Fig.  1 1)  is  more  sensitive  than  the  use  of  a  single  large  outer  window.  The 
streak  needs  to  be  detected  in  only  one  of  these  windows.  The  other  background  pixels 
do  not  contribute  to  diluting  the  streak  statistics.  With  this  filter,  the  logical  conditions 
to  detect  a  faint  star  are: 

condition  1 :  p;n  >  3  an  , 
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and 


condition  2:  max  (pout.  i )  <  2  an  , 

where  an  is  the  estimation  of  the  local  noise  standard  deviation,  |_i0ut-i  is  the  mean  signal 
evaluated  for  each  one  of  the  twelve  outer  windows  and  max  (|a0ut-i )  is  the  maximum 
value  of  these  estimates.  Note  that  the  threshold  values  ‘2an’  and  ‘3an’  have  been  set 
arbitrarily  and  could  be  fine  tuned  but  their  values  are  not  critical.  Also,  the  signal  has 
already  been  processed  (background  removed)  and  the  mean  background  value  is  not 
zero  because  the  pixel  values  are  now  only  positive  integers  with  values  starting  at  zero. 
The  idea  is  simply  to  define  a  criteria  like 

Pout-i  <  (Pn  +  tolerance  factor) 

where  pn  is  the  average  background  signal  where  non-zero  values  are  due  to  the  noise 
only.  The  use  of  an  is  very  practical  because  it  allows  the  algorithm  to  adapt  to  the 
image  noise  level.  Finally,  the  inner  threshold  value  has  to  be  higher  than  the  outer 
threshold  value  to  avoid  the  situation  where  a  satellite  streak  is  detected  in  the  inner 
window  and  not  in  the  outer  window;  like  a  star.  In  such  cases,  the  streak  would  be 
systematically  erased  pixel  per  pixel  by  the  algorithm.  An  inner  threshold  higher  than 
the  outer  threshold  forces  the  detection  of  the  streak  in  the  outer  windows  first. 

As  shown  in  Table  1,  this  filter  is  good  at  erasing  faint  stars.  It  is  also  good  at  preserving 
the  streak  pixels,  even  though  some  faint  stars  are  not  erased  to  ensure  preserving  that 
streak.  Usually,  those  undeleted  stars  are  faint  enough  to  create  no  significant  detection 
with  the  streak  detection  algorithm  described  in  Section  3.7. 

The  main  issue  occurs  when  the  bright  stars  have  larger  profiles.  Figure  3  shows  that 
very  bright  stars,  suffering  from  blooming  effects,  have  larger  profiles  that  may  overlap 
the  external  windows,  which  are  supposed  to  only  measure  the  background.  In  such 
cases  the  filter  will  fail.  However,  this  filter  can  be  adapted  to  have  two  or  even  three 
passes  with  different  parameters  for  window  size  and  threshold  values.  Hence,  for  bright 
stars,  the  inner  window  size  is  increased  to  1  lxl  1  or  even  13x13  pixels  and  the  threshold 
value  for  the  background  estimation  in  the  outer  window  is  also  increased.  Therefore, 
the  detection  conditions  become  something  like: 

condition  1 :  p;n »3on  (  ...  typically  >  20  an  ) 

(typically:  pm  >  1000  count,  i.e.  1%  to  2%  of  the  dynamic  range  of  the  CCD  camera)  and 

condition  2:  [  max  (pout.j  )  <  (p;n  /  10)  ]  or  [  max  (pout-i  )  <  2on  ]. 

Hence,  the  brighter  the  star  in  the  inner  window,  the  higher  the  outer  window  threshold 
is.  When  detecting  a  bright  star,  this  makes  the  filter  less  sensitive  to  background  noise. 
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Table  1.  Filter  behaviour  depending  on  the  cases 


Encountered 

situation 

Inner  window 

Pin 

Outer  windows 

max  (Pout-i ) 

result 

Faint  streak 

Pstreak  <  3(7n 

undetected 

Might  be 
detected 

No  action,  streak  is  not  erased 

bright  streak 

Pstreak  >  3(7n 

detected 

detected 

No  action,  streak  is  not  erased 

Very  faint  star 

Pin  <  3a  n 

undetected 

undetected 

No  action  but  such  a  faint  signal 
do  not  affect  the  streak  detection 
algorithm 

Faint  star 

Pin  >  3(7  n 
and  no  apparent 
blooming  neither 
apparent  decay 

detected 

undetected 

Requires  the  execution  of  the 
star-removal  processing 

Bright  star 

Pin  »  3a  n 
with  apparent 
decay  or  blooming 

detected 

detected 

No  action:  require  another  filter 
with  larger  inner  window.  The 
star  intensity  is  an  indication  of 
the  window  size  that  should  be 
used. 

Pair  of  faint  stars 

detected 

detected 

No  action:  this  is  the  rare  case 
where  the  filter  is  eluded. 

Star  close  to  a 
streak 

detected 

detected 

The  star  is  not  erased  and  the 
streak  remains  intact. 

Bright  stars  are  very  easy  to  detect  by  nature  and  the  threshold  values  are  not  critical. 
The  idea  remains  that  when  a  bright  object  is  detected  in  the  inner  window,  one  must 
verify  that  the  brightness  is  fainter  in  the  outer  window  to  avoid  detecting  and  erasing 
bright  streaks. 

In  summary,  the  diagram  of  Fig.  12  illustrates  the  entire  process  proposed  to  detect  and 
erase  stars.  It  begins  with  erasing  the  already  known  saturated  stars  that  were  detected 
before  the  removal  of  the  image  background  (  Section  3.4),  then  by  erasing  the  bright 
stars  and  finally  by  erasing  the  faint  stars.  Bright  stars  are  copied  and  preserved  in  a 
separated  image  for  telemetry.  Usually,  software  for  calculating  the  telemetry  works 
better  with  these  clean  images  where  all  the  artefacts  have  been  removed.  Faint  stars  are 
simply  reset  to  zero  because  they  are  completely  confined  in  the  inner  window.  Bright 
stars  with  profiles  of  various  widths  are  erased  using  a  special  process  explained  in  the 
next  section. 

7.2  Erasing  a  star: 

There  are  several  methods  for  erasing  a  star.  There  are  the  bad,  the  very  bad  and  the 
appropriate  method  to  do  it.  All  methods  are  based  on  the  same  trick,  i.e.,  when  a  star  is 
detected,  its  pixels  are  set  to  zero.  The  main  problem  is  to  determine  which  pixels 
belong  only  to  the  star,  which  ones  are  completely  outside  of  its  area  of  influence  and 


22 


DRDC  Valcartier  TR  2005-386 


which  ones  are  partially  influenced  by  the  star  and  that  may  also  contain  other 
information  (background  or  near  objects).  The  star  profile  is  unique  for  each  star.  It  is  a 
combination  of  optical  transfer  function  plus  the  degradation  produced  by  the 
atmospheric  transmission  and  turbulence  and  the  blooming  effect.  Then,  for  each  bright 
star,  the  profile  must  be  measured  within  the  image.  Various  approaches  or 
combinations  of  approaches  can  be  used;  Fig.  13  illustrates  an  example  where  the 
original  image  is  processed  by  two  different  methods  giving  different  results. 


Step  1 


Step  2 


Step  3 


Figure  12.  The  three  steps  of  the  star  detection  and  erase  process 
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A 


Unaffected  faint  star 


Zone  of  subtracted  profile 


Figure  13.  Star  removal  process:  A)  an  enlarged  view  of  Fig.  19C,  B)  image  where  the  bright  stars  have 
been  deleted  by  setting  their  bright  pixels  to  zero  and  C)  image  where  only  the  brightest  star  has  been 
erased  by  setting  their  central  disks  to  zero  and  by  subtracting  their  measured  profiles  over  a  larger  area. 

This  process  removes  the  crown  shape  artefact  that  remains  in  B 

The  issue  is  to  determine  which  pixels  must  be  erased.  The  first  method  consists  in 
setting  the  brightest  star  pixels  (above  a  certain  threshold)  to  zero  (as  illustrated  in  Fig. 

13B).  A  better  method  can  be  to  use  a  circular  mask  (of  a  predetermined  radius)  where 
all  the  pixels  contained  in  this  area  are  reset.  In  that  case,  the  center  of  the  star  is 
determined  by  its  pixel  of  maximum  intensity  and  the  circular  mask  is  centered  on  that 
pixel.  This  method  is  very  appropriate  for  faint  stars  where  only  a  small  window  needs 
to  be  erased.  This  method  is  used  in  step  3  of  the  diagram  in  Fig.  12. 

However,  this  method  is  not  appropriate  for  brighter  stars  (particularly  the  saturated 
ones,  such  as  those  considered  in  steps  1  and  2)  where  the  blooming  creates  larger  star 
profiles.  In  those  cases,  the  erasing  of  the  central  disk  leaves  a  crown-shape  artefact  (as 
illustrated  in  Fig.  13B).  This  is  not  desirable.  Larger  marks  can  be  used  for  brighter 
stars  but  this  solution  is  not  an  ideal  because  very  bright  stars  would  require  masks  so  big 
that  large  parts  of  the  image  would  be  erased  including  all  the  other  objects  and  perhaps 
the  sought  streak.  Therefore,  a  clever  method  is  required. 

For  erasing  bright  and  saturated  stars,  the  selected  method  combines  central  mask  and 
peripheral-crown  erasing  processes.  The  size  of  the  central  circular  mask  is  set  to  the 
size  of  the  detection  inner  window;  typically  a  diameter  of  13  pixels  (for  the  13x13 
window  of  the  step  2  of  the  diagram  of  Fig.  12).  This  size  can  be  increased  for  rare 
saturated  stars,  according  to  the  number  of  saturated  pixels,  which  is  an  indication  of  the 
importance  of  the  blooming  and  bleeding  effects.  However,  as  mentioned  previously, 
excessive  bleeding  indicates  that  the  image  is  too  corrupted  to  be  of  any  value  for 
automatic  processing. 

The  peripheral-crown  erasing  process  is  a  more  complicated  process  because  the  star 
profiles  p(r)  need  to  be  estimated  for  each  star.  Then,  this  profile  is  subtracted  from  the 
image  area  affected  by  the  star.  Here  is  the  process  in  detail. 
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The  central  coordinate  (k0,l0)  of  the  star  is  first  found  (the  brightest  pixel  if  sub-pixel 
accuracy  is  not  required).  For  all  pixels  Tkl’  in  the  area,  the  radial  distance  |  r  |  (in  pixel 
unit)  is  calculated  using: 

|  r  |  =  |  ?k!  -r0\  =  ^(k-kj2 +(l-lo)2  . 

Then,  for  all  pixels  ‘I’  of  a  similar  distance  ‘  |  r  |  ’  that  form  the  image  subgroup 

e{I(|  r  | )},  the  estimated  profile  <p(r)>  is  obtained  by  calculating  the  median  value  of 
that  group  of  pixels: 

<p{r)>  =  median  [  e  {I(|  r\  )  }  ]. 

Other  operators  were  tried  but  the  median  value  was  found  to  be  the  most  robust  in  the 
presence  of  noise  and  nearby  objects.  For  instance,  in  Fig.  13-C,  a  very  faint  star  is 
included  in  the  area  where  the  bright-star  profile  is  evaluated  and  subtracted.  As  one  can 
see  this  very  faint  star  does  not  affect  the  median. 

Finally,  this  profile  is  subtracted  from  the  image: 

I'ki  =Iki-<p(|  ra  ~r0\)>  , 

excluding  the  most  central  pixels,  which  are  systematically  set  to  zero  over  the  area  of 
the  circular  mask.  Also,  <p(  |  r  kI  —  ro  |  )>  requires  a  few  interpolations  of  the  recorded 
function  <p(r)>. 

The  result  is  an  image  where  the  bright  stars  disappear  without  affecting  the  nearby 
objects.  In  fact,  there  are  always  artefacts  generated  by  any  processing.  Fortunately  for 
this  case,  remaining  artefact  amplitude  is  comparable  with  the  background  noise  level. 

The  result  of  the  star  removal  process  can  be  seen  in  Fig.  19.  Figure  19-C  shows  the 
image  after  the  gradient  background  was  removed.  Figure  19-E  shows  the  same  image 
after  removal  of  the  saturated  and  brightest  stars.  Finally,  Fig.  19-G  shows  again  that 
same  image  after  removal  of  (Fig.  12,  step  3)  all  the  faint  stars.  In  that  last  case,  the 
image  contrast  has  been  boosted  and  one  can  see  that  only  the  faint  satellite  streak 
remains  with  a  few  pairs  of  closed  stars  that  succeeded  to  elude  the  double-gate  filters. 


7.3  Creation  a  clean  star  map 

One  of  the  important  steps  in  the  image  analysis  is  the  determination  of  the  image 
telemetry.  The  exact  positions  of  known  stars  and  pixel  equivalent  celestial  coordinates 
(in  right  ascension  and  declination)  need  to  be  evaluated.  This  is  done  with  commercial 
software  that  can  recognize  stars,  simply  by  using  rough  telescope  pointing 
measurements  and  a  star  database.  However,  this  software  frequently  fails  to  identify  the 
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stars  when  the  image  background  is  too  intense,  when  there  is  too  much  noise  or  simply 
when  there  are  too  many  stars  (in  this  last  case,  the  software  cannot  converge  toward  a 
unique  solution). 

The  performance  of  such  a  software  package  can  be  improved  with  a  clean  star  image. 
When  a  star  is  erased  from  the  main  image,  it  is  concurrently  copied  into  a  parallel  image 
to  save  telemetry  information.  By  doing  so,  one  generates  a  clean  image  without 
background  and  with  only  the  optimal  number  of  stars  (the  brightest  ones)  required  by 
the  telemetry  software.  The  image  of  Fig.  19-D  only  shows  the  brightest  stars  that  were 
copied  into  it.  This  image  was  processed  by  the  ‘CCDSoft’  software  (created  by  Bisk 
Software)  and  all  the  stars  were  recognized.  Their  magnitudes  are  listed  in  Fig.  19-F. 
With  this  information  available,  the  real  position  of  the  pixel,  in  terms  of  right  ascension 
and  declination,  is  known  with  precision  (about  one  arc  second  of  error).  The  recognized 
stars  provide  information  about  the  image  photometry,  which  can  be  useful  for  future 
analysis  of  the  satellite  streak  brightness. 
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8.  Streak  filtering  and  detection 


The  technique  for  detecting  the  satellite  streak  is  based  on  well-known  matched  filter 
techniques.  However,  these  methods  are  usually  very  sensitive  to  the  background  signal 
and  can  create  several  false  alarms.  This  is  why  the  previous  filtering  process  was  so 
critical  for  the  production  of  an  image  with  a  reduced  noise  and  bright  star  content. 

Now,  one  can  expect  good  performance  for  a  matched  filter  applied  to  such  an  enhanced 
image. 

8.1  The  a-priori  knowledge  and  the  design  of  the  matched 
filter 

To  design  a  matched  filter,  one  must  have  a  priori  knowledge  (template)  of  the  desired 
object.  This  is  exactly  the  case  for  the  current  problem.  Here  the  problem  is  not  to 
detect  any  kind  of  satellite  streaks  but  only  the  streak  of  the  satellite  that  must  be 
reacquired.  This  satellite  and  its  orbital  parameters  are  already  known.  There  is  only 
incertitude  in  its  exact  current  position  because  of  the  degeneration  of  its  orbit.  This  is 
why  it  must  be  periodically  re-acquired.  Hence,  at  the  moment  of  image  acquisition,  the 
satellite’s  position  can  be  predicted  (using  orbital  models)  at  the  time  at  which  the 
camera  shutter  is  opened  and  at  the  time  at  which  it  is  closed,  thus  providing  two 
positions  for  two  times  of  reference.  These  estimated  positions  should  be  close  to  the 
real  satellite  positions.  They  should  be  inside  the  telescope  field  of  view  (by  system 
design  and  operation  constraints).  With  this  information,  the  length  and  orientation  of 
the  satellite  streak  in  the  image  can  be  estimated.  The  only  remaining  unknowns  are  the 
real  position  and  the  brightness  of  the  satellite.  This  information  will  be  provided  by  the 
convolution  peak  (the  result  of  the  matched  filter). 

Figure  14  (in  the  left  frame)  shows  an  example  of  the  expected  satellite  streak.  It  is  in 
fact  modeled  by  the  line  segment  that  links  the  two  expected  satellite  positions  Pi  and  P2 
for  the  two  times  t  and  t+dt.  This  line  segment  almost  corresponds  to  the  matched  filter 
that  will  be  convolved  with  the  image.  The  matched  filter  must  not  inject  a  phase  shift  in 
the  solution  of  the  convolution.  Therefore,  the  filter  mask  must  be  correctly  generated, 
i.e.,  the  object  center  must  be  correctly  placed  on  the  window  origin  (i.e  the  zero 
coordinates)  such  as  indicated  in  Fig.  14.  In  other  words,  the  object  is  centered  on  the  top 
left  comer  of  the  window.  That  filter  is  also  normalized  to  generate  a  correlation  peak 
similar  in  intensity  to  the  object  found  in  the  original  image.  Thus,  the  amplitude  of  this 
generated  line  segment  is  1/n,  ‘n’  being  the  number  of  pixels  included  in  that  segment. 


8.2  The  matched  filter  technique 

Figure  15  illustrates  the  well-known  principles  of  the  matched  filter  technique.  The 
convolution  of  the  image  with  the  filter  can  be  performed  by  solving  the  classical  double 
integral  for  all  pixels.  This  technique  is  only  efficient  in  terms  of  processing  speed  for 
objects  of  limited  sizes.  For  long  streak  (often  more  than  one  hundred  pixels),  it  is  better 
to  use  the  global  process  based  on  Fourier  transform,  such  as  shown  in  Fig.  15. 
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Figure  14.  Representation  of  the  object  and  creation  of  the  filter.  In  the  filter ;  the  object  is  generated 
from  the  point  of  origin  in  order  not  to  induce  a  position  bias  in  the  response  of  the  convolution 


R:  Result  of  the 

I :  Input  image  S  :  Streak  to  detect  convolution 


I  =  Fourier  transform  S  =  Fourier  transform  R:  Result  of  the 

of  the  input  image  of  the  streak  filter  convolution 


Figure  15.  The  convolution  peak  is  obtained  by  first  calculating  the  Fourier  transform  for  both  the 
image  and  the  filter  (or  object),  by  multiplying  them  (with  the  complex  conjugate)  and  finally  by 
calculating  the  inverse  Fourier  transform  of  the  result. 

The  Fourier  filtering  technique  is  well  described  in  any  signal  processing  textbook.  For 
the  benefit  of  the  reader  not  familiar  with  that  technique,  here  is  a  brief  description.  Let 
us  call  the  image  acquired  by  the  CCD  and  pre-processed  to  remove  the  background  and 
stars  ‘I’,  its  Fourier  transform  (/  =  F  (I)  )  7  ’,  ‘S’  the  image  of  the  synthesized  streak,  ‘S’ 
its  Fourier  transform  (S  =  F  (S) )  and  ‘S  ’its  complex  conjugate.  The  convolution  of 
the  image  with  the  object  is  obtained  by  multiplying  their  Fourier  transforms  and  then  by 
calculating  the  inverse  Fourier  transform  of  the  result,  i.e.:  R  =  F _1  {  /  S').  The 
convolution  peak  that  appears  in  the  filtered  image  ‘R’  indicates  the  location  of  the 
searched  streak.  Figure  19-J  shows  an  example  of  such  a  convolution. 
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8.3  Analysis  of  the  shape  of  the  convolution  peak 


There  may  be  several  false  alarms  generated  by  bright  objects  that  do  not  really  match 
the  filter.  These  objects  have  so  much  energy  in  their  structures  that  even  a  partial  match 
is  still  more  intense  than  a  perfect  match  with  a  faint  streak.  This  is  the  case  of  Fig.  19-J 
where  three  pairs  of  relatively  bright  stars,  which  survived  to  the  star  removal  process, 
generate  three  convolution  peaks  more  intense  than  the  faint  satellite  streak.  Therefore, 
the  response  of  the  matched  filter  must  be  analyzed  further  to  see  how  the  real  match  can 
be  discriminated  from  the  false  alarms.  Fortunately,  the  answer  to  this  problem  is  quite 
simple  as  demonstrate  below.  Let  us  see  the  appearance  of  the  shape  of  the  convolution 
peak  for  both  objects  (streak  and  star). 

When  observing  Fig.  16,  one  can  see  that  there  is  a  major  difference  between  the  results 
of  the  matched-filter  for  a  streak  or  a  star.  In  theory,  the  convolved  streak  has  a 
pyramidal  shape  with  an  amplitude  ‘a’  (also  the  amplitude  of  the  streak  in  the  original 
image)  while  a  convolved  star  has  almost  a  rectangular  shape  (convolved  by  the  narrow 
star  Gaussian  profile)  with  an  amplitude  of  ‘2b/n’,  where  ‘b’  is  the  maximum  star 
intensity  and  ‘n’  the  streak  length  in  pixel  (or  the  length  of  the  line  segment  that  is  used 
as  template).  In  practice,  the  result  of  the  filtering  method  is  modulated  by  noise.  Also, 
the  pyramidal  shape  of  the  convolution  peak  may  be  distorted  by  the  satellite  fading 
(principally  due  to  satellite  rotation). 
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Figure  16.  Matched  filter  convolution  technique:  response  of  a  streak  and  a  star  when  convolved  by 

a  ‘Rectangle’  function 
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8.4  Iterative  matched-filter  technique 

This  observation  is  important  and  suggests  what  to  do  in  the  next  processing  step.  The 
shape  and  intensity  of  the  streak  convolution  peak  is  similar  to  the  original  streak.  This 
is  completely  different  for  a  star  where  the  convolution  peak  is  a  lot  fainter  than  the 
original  star. 

This  observation  suggests  the  following  processing  method.  For  each  detected 
convolution  peak,  the  ‘intensity’  information  can  be  used  to  clip  the  original  signal  and 
produce  a  new  enhanced  image.  In  this  new  image,  the  original  streak  is  preserved  while 
the  stars  are  attenuated. 

Here  is  the  basic  principle  of  the  clipping  process.  Each  object  detected,  indicated  by  a 
convolution  peak,  is  clipped  with  an  adapted  ‘clipping  mask’.  The  shape  of  this  mask  is 
the  streak  shape.  Its  intensity  is  that  of  the  convolution  peak.  Therefore,  the  streak 
clipping  mask  is  very  similar  to  the  real  streak.  On  the  other  hand,  the  star  clipping 
masks  are  larger  and  fainter  than  their  corresponding  stars.  The  clipping  process  consists 
in  calculating  the  minimum  between  the  original  image  and  the  clipping  masks.  Hence, 
the  streak  remains  unchanged  while  the  stars  are  severely  truncated.  Such  an  enhanced 
image  is  shown  in  Fig.  191. 

Since  the  streak  remains  unchanged  after  this  process,  this  new  enhanced  image  can  be 
processed  again  with  the  matched  filter.  This  will  produce  a  second  generation  of 
convolution  peaks  where  the  streak  convolution  peak  is  unchanged  but  where  the  star 
convolution  peak  is  fainter  (Fig.  19  L).  This  new  result  can  be  used  to  generate  the 
second  generation  clipping  masks  and  the  entire  process  repeated  again  and  again. 

Figure  19J  and  19L  show  the  convolution  peaks  after  one  and  three  iterations  (with  the 
other  method  described  below)  while  Fig.  191  and  19K  show  the  clip  images  after  as 
many  iterations.  Experiments  (Ref.  10)  showed  that  12  iterations  produce  the  best 
results.  There  is  no  visible  improvement  with  additional  iterations. 

This  technique  appears  very  difficult  to  implement.  However,  there  is  a  more  practical 
method  to  achieve  the  same  result.  Rather  than  searching  for  every  convolution  peak 
and  generating  several  intensity  masks,  one  can  simply  directly  use  the  previous 
convolution  result.  The  convolution  peaks  of  the  stars  already  have  the  appropriate 
shape  (almost  a  rectangle  function  such  as  shown  in  Figures  16  and  17),  location  and 
intensity.  Therefore,  the  convolved  image  can  be  used  as  a  global  intensity  mask,  which 
includes  all  objects  at  the  same  time.  This  is  good  for  the  stars  but  not  for  the  streak 
since  this  is  the  only  object  with  a  triangular-shape  convolution  peak.  However,  by 
multiplying  the  convolved-image  intensity  by  two,  such  as  shown  in  Fig.  17,  this 
triangular  peak  completely  contains  the  streak  profile.  Figure  21  shows  the  same  thing 
but  with  real  data.  This  convolved  image  (multiplied  by  2)  can  be  used  as  a  clipping 
mask.  With  this  simpler  clipping  method,  the  stars  are  less  severely  clipped,  whereas  the 
streak  is  entirely  preserved. 

Figure  17  illustrates  the  original  signal,  the  shape  of  the  convolution  peaks  (2x)  for  the 
streak  and  for  a  star,  along  with  the  resulting  clipped  signal.  Figure  18  illustrates  the 
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iterative  process  that  involved  all  the  Fourier  transform  and  signal  clipping.  As 
mentioned  above,  the  whole  process  practically  converges  after  three  iterations.  After 
that,  the  iterative  matched  filter  technique  provides  no  significant  improvement. 
Flowever,  an  ongoing  study  (ref.  10)  shows  that  the  best  SNR  is  reached  after  12 
iterations. 


Convolution  peak 
multiplied  by  2 

Streak  profile 


-n  n 


Minimum  of  both  functions: 


—  a 

1 

-n/2 

The  streak  profile  remains  unchanged 


Minimum  of  both  functions: 

-4b/n 


-n/2  11/2 

The  star  profile  is  severely  truncated 


Figure  17.  Clipping  the  original  image  using  a  local  intensity  threshold  function ,  which  is  made  using 
the  result  convolution  of  the  original  image  with  the  matched  filter. 
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Figure  18.  Iterative  matched-filter  technique 
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9.  Result 


The  results  of  all  the  processing  methods  explained  above  are  illustrated  in  Fig.  FI  9. 
Starting  with  the  raw  image  (Fig.  19A),  the  background  is  removed  (Fig.l9C),  the  stars 
are  erased  (Fig.  19G)  and  the  satellite  streak  is  extracted  (Fig.l9L).  The  images  speak  by 
themselves;  this  processing  is  extremely  efficient. 

The  processing  methods  are  very  efficient  in  terms  of  computer  requirements.  The  first 
implementation  in  Matlab  was  sluggish.  Flow  ever,  after  once  recoded  in  C++,  the 
complete  program  was  executed  in  only  a  few  seconds  with  an  old  processor  (500  MFlz). 
It  is  even  faster  to  process  images  than  to  download  them  from  the  CCD  camera. 
Therefore,  in  this  case,  the  processing  time  is  not  a  limitation. 

Also,  the  information  of  the  satellite  streak  is  well  preserved.  Figure  20  shows  an 
enlarged  view  of  the  satellite  streak  1)  after  the  removal  of  the  background  and  2)  at  the 
end  of  the  process.  Stars  close  to  streak  were  erased  and  the  noise  was  almost 
completely  eliminated  (particularly  by  the  clipping  process  during  the  steps  of  the 
iterative  matched  filtering).  Even  the  star  superposed  on  the  streak  was  attenuated  by  the 
clipping  process,  which  improves  the  perception  of  the  streak.  However,  the  streak  was 
shortened  by  a  few  pixels,  which  indicates  that  the  clipping  process  (in  the  iterative 
matched  filter)  could  be  improved.  Since  the  streak  signal  is  close  to  the  noise  level 
(SNR=2  in  the  present  case),  the  convolution  peaks  are  noisy  and  so  is  the  clipping 
mask.  However,  the  profiles  showing  these  thresholds  in  Fig.  21  seem  to  indicate  that 
the  methods  perform  as  predicted  in  Fig.  17.  The  method  also  seems  to  be  robust  to 
noise. 

The  methods  presented  here  were  tested  with  success  on  several  images.  Not  every 
possible  situation  has  yet  been  tested,  but  up  to  now,  detections  are  not  problematic  and 
the  algorithm  consistently  provides  good  performance.  For  one  of  the  worst  cases  where 
the  SNR  is  only  0.75  (Fig.  4C)  and  the  noise  standard  deviation  is  only  3  counts  per  pixel 
(i.e.,  the  streak  is  only  2  or  3  counts  above  background)  the  streak  was  clearly  detected. 
The  problem  with  such  a  signal  is  the  streak  analysis  afterwards.  Even  though  it  is 
detected,  its  telemetry  (positions  of  the  end-points,  center  of  mass,  etc.)  is  very  uncertain. 
However,  the  processing  for  such  analysis  requires  more  elaborate  numerical  analysis 
methods,  which  was  not  the  goal  of  the  present  study. 
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Figure  20.  Satellite  streak  before  the  matched  filter  (but  with  background  removed  and  where  the 
brightest  star  have  also  been  removed)  and  that  same  streak  after  noise  smoothing  and  three 

matched-filter  iterations. 


Figure  21.  Profiles  extracted  from  Figs.  19  G  and  J  showing  the  satellite  streak  and  a  star  (that 
survives  to  the  star  erasing  process)  before  and  after  convolution.  This  figure  (like  Fig.  17)  illustrates 
with  real  data  how  the  clipping  (achieved  through  the  convolution  of  the  image  with  the  matched 
filter)  preserves  the  streak  but  severely  truncates  the  star  pixel  value. 
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10.  Future  works 


The  algorithm  described  previously  performs  very  well.  However,  there  is  still  room  for 
improvement.  For  example,  some  parameters  could  be  adjusted  to  make  the  algorithm 
more  robust  in  the  presence  of  noise.  For  instance,  the  threshold  mask  could  be 
redesigned  to  better  preserve  the  streak  pixels.  For  example,  a  bias  of  an  could  be  added 
to  the  threshold  technique.  This  could  help  to  preserve  the  detected  streak  and  would 
leave  the  efficiency  of  the  filter  virtually  unaffected  when  a  star  is  to  be  truncated  (as  in 
Fig.  21).  However,  this  requires  more  experiments  that  will  be  made  in  the  next  phase  of 
development. 

It  would  be  interesting  to  add  post-detection  analysis  algorithms  in  order  to  reject  false 
alarms,  validate  the  detection,  extract  useful  information  and  generate  automatic 
detection  reports.  For  very  faint  streaks,  the  presence  of  false  alarms  becomes  an  issue. 
Very  often,  it  is  obvious  that  the  brightest  detection  peak  is  created  by  one  of  the 
remaining  bright  stars,  or  more  probably  by  a  pair  of  faint  stars  that  eluded  the  star 
erasing  process.  In  those  cases,  these  false  alarms  could  be  rejected  by  performing  basic 
morphological  analysis.  For  example,  the  shape  of  a  group  of  pixels  (indicated  by  the 
brightest  detection  peak)  could  be  evaluated  and  if  the  shape  does  not  correspond  to  the 
expected  streak  shape  (as  a  compact  group  of  pixels  created  by  a  star),  then  this  alarm 
can  be  discarded  and  the  analysis  could  continue  with  the  next  detected  object. 
Experiments  have  shown  that  this  kind  of  situation  occurs  when  the  streak  SNR  is  below 
2an,  where  the  valid  detection  is  very  often  the  second  or  third  detection  peak  in 
intensity.  Another  interesting  post-detection  analysis  method  is  the  extraction  of 
information  about  the  detected  streak.  Its  length,  orientation,  position  of  its  geometric 
center,  positions  of  its  end-points,  its  magnitude,  etc.  are  all  information  that  can  validate 
the  detection  and  allow  declaring  that  the  algorithm  succeeded  or  failed.  Subsequently, 
detection  reports  to  this  effect  can  be  generated. 


Note:  At  the  time,  this  report  is  being  published,  the  proposed  works  in  the  previous 
paragraph  were  undertaken  and  the  results  will  be  published  in  Ref.  10. 
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11.  Conclusion 


The  cascade  of  algorithms  described  in  this  document  were  implemented  and  tested  in 
Matlab  and  have  proven  to  perform  very  well.  A  more  optimized  version  was  developed 
in  C++,  which  has  also  proven  to  be  very  efficient  from  the  point  of  view  of  CPU  time. 
The  CPU  time,  required  to  process  an  image  for  satellite  streak  detection,  is  less  than  1 0 
seconds.  If  the  sensor  is  expected  to  acquire  a  new  observation  every  three  minutes  (this 
is  expected  for  the  CASTOR),  this  is  less  than  the  spare  CPU  time.  Therefore,  the 
processing  time  is  not  an  issue  for  real  time  operations.  The  images  can  be  processed 
immediately  after  their  acquisition. 

The  algorithms  were  implemented  in  the  automatic  processing  software  (at  the  SOC  in 
DRDC  Ottawa)  for  the  operation  of  the  three  CASTOR  stations.  Up  to  now  (after  several 
months  of  operation  and  applied  to  hundreds  of  images),  it  has  proved  to  be  very 
efficient  and  seems  to  be  able  to  detect  everything  with  a  SNR  above  2  with  almost  no 
false  alarms.  However,  due  to  the  number  of  pixels  in  the  streak  (hundreds  of  pixels 
with  SNR>2,  i.e.,  the  signal  is  still  strong),  it  is  expected  that  the  method  can  still  be 
improved  to  detect  fainter  streaks.  For  fainter  streak  intensity,  the  actual  algorithm 
performance  is  limited  by  false  alarms  (not  by  the  sensitivity)  and  the  suggestion  made 
above  (about  morphological  analysis)  may  be  a  solution  to  lower  the  detection  limits. 

Up  to  now,  detections  have  been  achieved  with  a  pixel  SNR  lower  than  one.  At  this 
level,  it  is  however  difficult  to  extract  additional  information  (like  the  end-point 
coordinates)  about  the  streak.  Historically,  the  acquisition  and  analysis  are  performed 
with  a  human  in  the  loop.  One  single  shot  is  acquired  and  the  analyst  tries  to 
characterize  the  single  streak  by  measuring  its  end-points.  For  very  low  SNR,  end-point 
cannot  be  determined  with  precision  and  the  analyst  can  only  report  the  observation  with 
inaccurate  details.  By  using  numerical  analysis,  it  is  possible  to  extract  more  reliable 
information  such  as  the  measurement  of  the  geometric  center  of  the  streak,  which  is  less 
sensitive  to  the  noise  than  the  end-point  estimation.  Furthermore,  an  automatic  system 
could  acquire  a  strobe  streak  (with  multiple  frame  expositions),  providing  several  streak 
segments,  which  would  allow  one  to  calculate  several  geometric  centers  and  verify 
whether  these  detections  are  coherent.  Hence,  the  processing  could  increase  the 
sensitivity  and  reliability  of  the  detection  system,  but  it  would  require  changing  the  data 
acquisition  procedure.  The  detection  algorithm  presented  in  this  document  was  proven 
to  perform  very  well.  This  tool  is  now  being  considered  for  inclusion  within  the 
operation  tools  of  the  observation  stations. 
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List  of 

symbols/abbreviations/acronyms/initialisms 


DND 

Department  of  National  Defence 

SofS 

Surveillance  of  Space 

soc 

Sensor  Operation  Control 

ssoc 

Surveillance  of  Space  Operation  Center 

PSF 

Point  Spread  Function,  also  called  the  optical  impulse  response 

SSN 

Space  Surveillance  Network 

CASTOR 

Former  name  of  the  ground  stations 

SofS/CD 

SofS  Concept  Demonstrator;  new  name  of  the  ground  stations 

CCD 

Charge  Couple  Device 

SNR 

Signal-to-Noise  Ratio 

ADC 

Analog  to  Digital  Converter 
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